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Abstract 

Eye movement data are outputs of an analyser tracking the gaze when a person 
is inspecting a scene. These kind of data are of increasing importance in scientific 
research as well as in applications, e.g. in marketing and man-machine interface 
planning. Thus the new areas of application call for advanced analysis tools. Our 
research objective is to suggest statistical modelling of eye movement sequences 
using sequential spatial point processes, which decomposes the variation in data 
into structural components having interpretation. 

We consider three elements of an eye movement sequence: heterogeneity of the 
target space, contextuality between subsequent movements, and time-dependent 
behaviour describing self-interaction. We propose two model constructions. One 
is based on the history-dependent rejection of transitions in a random walk and 
the other makes use of a history-adapted kernel function penalized by user-defined 
geometric model characteristics. Both models are inhomogeneous self-interacting 
random walks. Statistical inference based on the likelihood is suggested, some ex¬ 
periments are carried out, and the models are used for determining the uncertainty 
of important data summaries for eye movement data. 
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1 Introduction 


Eye movements reflect brain functions, revealing information on ongoing cognitive pro¬ 
cesses, and can be recorded by eye trackers in a cost-efficient way. Eye movement data 
are spatio-temporal and consist of time sequences of fixations, points in the target space 
where the gaze stays for a while, and of saccades, which are rapid jumps between hxations. 
An example of eye movement data can be seen in Figure 



Figure 1: First 100 hxation points of one subject on a painting called Black Bow (1912) 
by Wassily Kadinsky. The arrows show the movement of the gaze during the hrst three 
seconds. 

Fixation locations in eye movement data are point patterns which can be modelled by 
means of spatial point processes. Point process statistics is a well-developed branch of 


spatial statistics increasingly used in applied sciences, see e.g. Illian, Penttinen, Stoyan, 


and Stoyan (2008) and Diggle (2013). Extensive software spatstat (Baddeley, Rubak, 


and Turner, 2015) has made efficient point pattern data analysis attractive. Point process 


statistics has been applied for eye movement data by (Barthelme, Trukenbrod, Engbert)] 
and Wichmann (2013), who use the spatial inhomogeneous Poisson point process to 


predict the hxation locations. The approach by Barthelme et ah (2013) aggregates the eye 


movement data over time but omitting all dynamics. Engbert, Trukenbrod, Barthelme, 


and Wichmann (2015) present a dynamical model that takes spatial interaction into 


account, but their model validation is based on characteristics of spatial point processes. 
A step towards a dynamic model is to add the temporal order of hxations, which leads 
us to the class of (hnite) sequential spatial point processes, see van Lieshout ( 2006a||b 
2009). If, in addition to the order, the time instances of occurrences of the points are 


recorded and included in the model, then the underlying process is a spatio-temporal 


point process, see 

e.g. 

Diggle, Kaimi, and Abellana 

(2010 

), and an application to eye 

movement data by 

Ylitalo, Sarkka, and Guttorp 

(^201f 

)). 



We consider eye movement data to be a realisation of a sequential spatial point process 


which allows us to extend the approach by Barthelme et ah (2013) for detecting new im¬ 


portant dynamic structures of data. The advantage of this approach is that the likelihood 
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is tractable and the simulation of realisations is straightforward. In addition, sequential 
point process modelling is a construction step for spatio-temporal point processes. 


For eye movement data, three structural components of sequential spatial point processes 
are central. Spatial heterogeneity of hxation pattern means that some parts of the scene 
get the observer’s attention more than others. This strong component is present in almost 
all eye movement data. It is usually modelled through a saliency map, which is calculated 
from the features of the scene (see e.g. Itti and Koch, 2000, Kiimmerer, Wallis, and Bethge 


2014), such that the most salient areas are expected to obtain more hxations. Dynamic 


contextuality is a saccadic property which describes the (metric) length of a jump from the 
current fixation to the next one: for example, nearby sites may be more favourable than 


the more distant ones (see e.g. Tatler, Baddeley, and Vincent, 2006). Both heterogeneity 


and contextuality are well-established in eye movement studies (see e.g. Barthelme et ah 
2013, Engbert et ah, 2015, Kiimmerer et ah, 2014). 


However, our empirical evidence shows that these two components are not sufficient, 
since e.g. they cannot model the learning effect. Similar hndings are made by |Engbert| 
et ah (2015) and Kiimmerer et ah (2014). Thus we are looking for simple mechanisms 


which could utilize the long-term dependence indicating the learning process during an 
experiment. One such possibility is self-interaction, which modihes the individual moves 
(saccades) by means of the history of the eye movement sequence. As an illustration, 
the observer prefers to inspect the whole scene at the beginning of the experiment and 
gradually focuses on a few details (see e.g. Locher, Gray, and Nodine, 1996). Here, 


we offer a tool for studying the self-interaction effect in eye movement sequences. We 
suggest new models for eye movement data to deduce the effect of structural components 
and to evaluate statistical variation in problem-specihc functional summary statistics. 
The suggested models have potential use also beyond the eye movement research, e.g. in 
ecology for modelling animal movements, and in user-interface studies. 

Our idea is that the history of the sequence changes the dynamics during the evolution 
of the eye movement sequence. We present two general principles for model construction, 
both of which are generalisations of the random walk in heterogeneous media. The hrst 
principle is the history-dependent thinning of transitions, which assigns smaller weights 
for suggested transitions being at odds with the chosen functional summary characteristic 
conditional on all previous hxations. This type of penalization is similar to the area- 
interaction process by Baddeley and van Lieshout (1995) in point process statistics. The 


second principle resembles the ARCH (autoregressive conditional heterogeneity) model, 
commonly applied in econometric time series analysis (Engle, 1982). Based on these 


principles, we present the construction of the two processes, how to simulate them, and 
how to estimate the parameters by the maximum likelihood method. Several summary 
statistics, assisted by Monte Carlo simulation, are applied in model evaluation. 


The new models are mainly of “statistical” nature, which means that they do not mimic 
the neural process, but they can capture essential variation in eye movement data. We 
will not employ all the generality the suggested new models are able to achieve. Instead, 
the objective is to present the new ideas in terms of rather simple models which still are 
useful in eye movement data analysis, especially in the study of the learning mechanism 
during an experiment, and in the derivation of statistical variation of important data 
summaries. Furthermore, this new approach will bring eye movement data analysis closer 
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to statistical inference. The motivation is the complexity of eye movement data, which 
are inhomogeneous in space and time, and the use of asymptotic inference, for example, 
is difficult to justify. 

The paper is organized in the following way. In Section 2, two new sequential spatial 
point process models are suggested. Simulation algorithms are given in Section 3 with 
simulation experiments demonstrating the models. In Section 4, eye movement data in 
the held of art study is modelled by using the new approach to deduce self-interaction. 
Section 5 contains some concluding remarks. 


2 Finite sequential spatial point process models 

Suppose l^n = (a^i,... ,Xn) is a sequence of time-ordered points in a bounded window 
IT C The corresponding unordered point set {xi,... is denoted by If 

(IT"', W") stands for the n-dimensional space of ordered points provided with the Borel a- 
algebra in IT", the density function / w.r.t. the Lebesgue measure is dehned sequentially 
as follows: fi{xi) stands for the probability density of the hrst point, and the conditional 
density of a further point x given k = (a^i, • • •, Xk), fc = 1,..., (n — 1), is denoted by 
fk+i{x\~^k), a transition probability density for the transition k ^ , x), and the 

joint density of is 


n—1 

fC^n) = fl{Xl) n fk+l{Xk+l\ltk). (1) 

k=l 

The density fi{x) can be assumed to follow a function (e.g. the saliency map) exposing 
the focal areas of the target, or modelling can be conditional on the hrst observation xi. 
The transition densities fk+i{x\l^k)-, k = 1,..., n — 1, rehect saccadic features of the eye 
movement sequence. 

A simple model for the transition would be a random walk in heterogeneous media, which 
is dehned as 


fk+iix^^k) oc a{x) K{xk,x), (2) 

where a{x) is non-negative and bounded in IT, and K{xk,x) is a Markovian kernel, i. e., 

K{xk,x) > 0 for all Xk,x G IT, and 
/ K{xk, u) du = 1 for all G IT , 


w 


k = 1,... ,n — 1. In this simple model, a{x) describes heterogeneity of the scene. It 
can be a known saliency map, an empirical saliency map estimated as the intensity of 


repeated hxation patterns (see the discussion in Diggle, Gomez-Rubio, Brown, Chetwynd, 


and Gooding! 2007|), or a model based prediction of the saliency map, for example, 

p 


a{x) = h I ) > 
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where the variables Zj{x) are the values of p feature vectors at x extracted from the scene 
using machine learning techniques, hj-.s are regression coefficients, and h is an adequate 
non-negative function (see e.g. Barthelme et ah, 2013). The Markovian kernel K{xk,x) 


describes the contextuality of subsequent fixations related to jump lengths. An example 
is the truncated Gaussian kernel 


K{xk,x) oc e , Xk,xeW, 


(3) 


where | Ixfc —a;| | is the Euclidean distance between the points Xk and x. For the rectangular 
window W = [a, b] x [c, d] the normalization term for ([^ can be written as 

27ra^{^{b) - - $(c)). 


where <1> is the c.d.f. of the standard normal distribution. This kernel penalizes large 
jumps and keeps the process inside the specified window W. Thus the model ([^ captures 
heterogeneity in the target and models the transitions (or saccades) in a Markovian way. 


The transition mechanism ([^ may be insufficient if data contain learning. For instance. 


a two-stage model describing the nature of an aesthetic experience (see e.g. Tocher, 2006 


Tocher et ah, 1996, 2007) suggests that a picture is first inspected globally and, after hav¬ 


ing obtained a gist of the scene, the viewer starts to concentrate on some details. Also, the 
visual information gathered from the scene affects our cognitive processes and attention, 
which again affect the movement of the gaze. By keeping these complexities of human 
attention in mind, we develop two models which try to catch the sequential adaptation 
in the eye movement sequence in a tractable manner using geometric reasoning. 


2.1 Self-interaction due to history-dependent rejection model 

First, we define a history-dependent rejection model (later: rejection model), in which the 
self-interaction mechanism is created by a reweighting probability function. This model 
penalizes the location of the next point in terms of coverage or recurrence composed by 
the previous points: the density for the transition k ^ (^fc , x) is assumed to be 

fk+i{x[^k) oc a(x) K(xk,x)7r(x,S(~^k,x)), (4) 

where 7r(x, S{ltk,x)) is the reweighting probability of x when proposed according to the 
density proportional to a{x)K{xk, x). Here, S{l^k,x) = S{xi,... ,Xk,x) is a measure of 
coverage or recurrence of the ordered sequence (xi,... ,Xk,x). 

Reasonable choices of the reweighting probability vr in the eye movement context are 
given below. 

Coverage-based reweighting 

Two measures of the coverage of a point set are the area of its convex hull and the area 
of the associated ball union. From now on we assume that the scene W is convex. The 
convex hull of a point set {itk}, denoted by ConvC^^), is the minimal convex subset of 
W which contains all the points of {"affc}- The convex hull is unique and invariant under 
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permutation of the points; hence it is the same for ordered and unordered sets. Let us 
denote by Sc{~^k) the area of ConvC^fc) and call it convex hull coverage. 

The ball union measure of a point set {1^k} is defined as 

k 

Bcov(l^fc) = \Jh{xi,r) n W, 

i=l 

where b{x,r) stands for the ball with radius r and centred at x. It is a “regionalized” 
version of the point set, where the r close neighbourhood of a point is taken into account, 
and which again is invariant under permutation. Its area SsC^k) is called ball union 
coverage. 

The rationale behind model (|^ is that the kernel function generates random jumps 
and the reweighting probability determines which of the proposed jumps are accepted, 
depending on the current coverage of the sequence and on the new suggestion. Con¬ 
sider the convex hull coverage first: if the new suggestion x is not in ConvC^fc), the 
odds ratio of acceptance w.r.t. a proposal y G Conv(”^fc) with ||x — XfcH = ||r/ — is 
Sc{'^k,x)/Sc{'^k)- 


A reasonable and simple choice of the geometric nature for the reweighting probability 
would be 


7r{x,S{'^k,x)) 


1 if a; e IT \ Conv( affc) 

P if X G ConvCaffc) 


(5) 


where p G [0,1] and k>l. If p = 1, we have the random walk model. When p < 1, this 
choice encourages locations outside the convex hull of previous points leading to faster 
coverage. When the convex hull of the points covers almost the whole scene, the process 
behaves like a random walk. The density for the transition itk , x) with the 

truncated Gaussian kernel is 


/fc+i(x|^fc) oc a(a;)e """" (lw\Conv(atfc)(a^) + Plconv(atfe)(a^)), (6) 

where P = 1, 2,..., n — 1, and !(■) is the indicator function. 


The convex hull coverage can be replaced by the ball union coverage. It is not as sensitive 
to distant points as the convex hull coverage but reacts to the “holes” in the point pattern. 
If Conv(^fc) is replaced by BcovC^fc) in the reweighting probability (|^, the process 
favours locations away from the previous points and hence reduces clustering if p is 
small. It should be noted that the ball union coverage measure requires the radius of the 
ball. 


Recurrence-based reweighting 

As a measure of recurrence we propose the number of earlier visits in a ball h{x, r) around 
a site X, formally SR{~^k, x) = '^b(xi,r){.x). However, instead of using all the previous 

points, at step k the number of earlier visits is calculated from the point set {'^k-i} 
omitting the most recent point Xk- This delayed recurrence measure 

fe-i 

x) = '^ lb{xi,r)ix) (7) 

i=l 


6 


is less confounded with the Markovian kernel K{xk,x) than SR{l^k,x) and is therefore 
used in this paper from now on. Note also that the recurrence measure is not invariant 
under random permutation. 


A simple model for the reweighting probability is 


7l{x, SR(l^k,x)) 


6 if SRiltk, > 1 
1-9 ifSR{^k,x)=0 


( 8 ) 


where 6 G [0,1] and k > 2. The odds ratio for accepting a location close to the points 
of {1^k-i} against accepting a location from an empty area is 0/(1 — 9). If 9 is close to 
1, the process favours clustering, and if 9 is small, the process avoids previously visited 
local areas around the points If 0 = 0.5, the process is a random walk. The 

density for the transition k ^ , x) with the truncated Gaussian kernel is 

fk+i{x\^k) oc a{x) e~^ + 0 (9) 

where A; = 2, 3,..., n — 1, (and f 2 {x\l^i) = fi{x)). 


In particular, the model defined through (|^ is among the simplest ones which satisfy 
our requirements of self-interacting nature. Note that we need the normalized transition 
kernel in the likelihood, because the scaling factor contains the parameters of the model. 
The normalizing integral can be computed using numerical integration. Its evaluation 
can be avoided in the simulation of the process, however. 


2.2 Self-interaction due to history-adapted model 


The motivation behind the history-adapted model arises from the saliency map idea and 


the two-stage model by Tocher and colleagues (Tocher, 2006, Tocher et ah, 1996, 2007). 


Heterogeneity of the target plays the main role at an early stage of the process evolution: 
the areas with high saliency (or intensity) get more fixations than the low-saliency areas. 
However, when the target has been inspected well enough, the jump lengths get shorter 
as if the process were mimicking a local inspection process. 


This model construction is intended for coverage type self-interaction. We apply directly 
an adaptive Markovian kernel K^i^{xk,x), where cj)k is a function of the points and 

determines the width of the kernel. Hence, the kernel changes in time and affects the 
jump lengths. The transition probability density can be written as 


fk+i{x\'^k) 

4>k 


a{x) K^^{xk,x) 
Jwa{u)K^^{xk,u)du 
(pki'^k) oc H{S{ltk)) 


( 10 ) 

( 11 ) 


where H{s) is decreasing in s. Here a{x) controls the heterogeneity of the target as 
in the rejection model, whilst H{s) models the progress of the coverage. This model 
resembles the autoregressive conditional heterogeneity model (ARCH) commonly applied 
in time series analysis for modelling volatility (Engle, 1982). While ARCH models use the 
information from q lagged values, our model is allowed to use the entire history. Again, 
this model is a self-interacting random walk. 
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In what follows, we make use of the specihc model 


fc) 


cx e ^ 

^ ^^-KSil^k)/\W\ 


( 12 ) 

(13) 


T,K > 0 and Xk,x G W, where S{ltk) is the coverage of If k = 0, the process 

is a random walk since the kernel does not change in time. This model contains two 
parameters, r describing the initial kernel width and k modelling the decay as a function 
of coverage. The transition is determined by the conditional density (10). Both convex 
hull and ball union coverages are suitable for this construction. 


2.3 Model fitting and statistical inference 


2.3.1 Model fitting 


We assume that an ordered sequence = {xi,, Xn) is observed in W. In what follows 
we suggest parameter estimation for the two models defined by (rejection model), and 


by (12) and (13) (history-adapted model) assuming that the non-negative heterogeneity 
component a{x) is fixed. In practice, the estimation of a{x) is problematic, as we have 
pointed out in Discussion. 

History-dependent rejection model 

The log-likelihood for the general rejection model (|4) is now a function of the model 
parameters. For the two parameter model defined by (^ the expression 


n—1 


1 


n—1 


= ^^og{a{xk+i)) - ^^\\xk-Xk+i 


k=l 


(14) 


k=l 


n—1 


log(p) X^lconv(ltfe)(a;fc+l) 


k=l 


n—1 „ 

- ^log f 

^=1 w 


a[u)e 


\Xk-u\\ 


(lw\Conv(^fc)('w) + plconv(^fc)('w)) dw 


is obtained. Here we use the convex hull coverage in the reweighting probability, but 
also the ball union coverage could be used. The log-likelihood function for the rejection 
model with recurrence ([^ is shown in Section 4.1, formula (16). The logarithm of the 
normalizing factor (the last line of (|I4))) can be computed by numerical integration. 
The optimization of l{a‘^,p) w.r.t. can be conducted using numerical optimization, or 
alternatively, one can solve the exponential family likelihood equation 


n—1 


n—1 


[\\xk - l^fc) = '^\\xk - Xk+i\ 


k=l 


k=l 


where the expectation is over the conditional distribution of a new random point U 
from the distribution fk+i{x\l^k)- This can be computed using Monte Carlo maximum 
likelihood (MCMCML, see [Geyer ( 1991[ )). 
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Maximizing the log-likelihood ( |I4| ) (or (16)) is costly due to the normalizing integral, 
and we have chosen to use the profile likelihood approach (see e.g. Davison, 2008, p. 
127): First, cr^ is solved for fixed p resulting in cr^(p). Then, cr^ in the log-likelihood is 
substituted by a‘^{p) giving the log-profile likelihood lp{p) = l{a‘^{p),p), which is then 
maximized w.r.t. p and p is obtained. Finally, we fix p = p and compute the correspond¬ 
ing value of (j^, namely = cr^(p). These steps of coordinate descent are iterated to 
obtain the maximum likelihood estimates. Alternatively, numerical optimization methods 
can be used. 


History-adapted model 

The kernel width 0^ of the fcth transition of the random walk is a function of the model 
parameters and is adapted to the history of the sequence through the coverage measure 
The log-likelihood for the general model given by (10) and (|TT|) is 


n—1 

k=l 


log{a{xk+i))+ logK^^{xk,Xk+i)-log J a{u)K^^{xk,u)du 

w 


In the special case of (12) and (13) the log-likelihood is 

n—1 n—1 


^(d«^) = X]log(«(3^fc+i)) - X] 


k=2 

n—1 


k=2 


20(r, k) 


\^k I 


10 - ± n 

— ^log / a{u)e~^ 

^=2 \Y 


,k) 


\\xk-u\\ 


du 


(15) 


with t) = re k = 1,... ,n — 1. 

The log-likelihood can again be maximized directly, or alternatively, the likelihood equa¬ 
tions can be derived and solved: the estimation equations are 


n—1 n—1 

'^Er^^dixk - u\\^\itk)/4>k = iixfc - xk+iW^/ipk 

k=2 k=2 

n—1 n—1 

{S{ltk)\\Xk-U\\^\ltk) /<Pk = E S'(^fe) Wxk-Xk+iW^/ipk, 

k=2 k=2 

where the expectations are over the conditional distribution fk+i{x\l^k) with parameters 
r and k. The estimation equations are in accordance with the maximum likelihood 
equations for the exponential family of distributions. 


2.3.2 Model evaluation 

Model evaluation of spatial dynamic models is typically based on the selected functional 
summary statistics which measure different features of the model, such as coverage, re¬ 
currence and jump length as a function of time (or order). In addition, a saliency map 
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is plotted together with the fixation locations. The random variation of the summary 
statistics is estimated from simulations. 


Model evaluation is done by estimating several summary statistics from data and plot¬ 
ting the estimates together with the model based simulated pointwise envelopes being 
a parametric bootstrap method (see e.g. Efron and Tibshirani, 1994, p. 53). These 


envelopes indicate statistical variation in the summary statistic under the parametric 
model assumption. It should, however, be noted that when using the pointwise envelopes 
as statistical tests, the multiple testing problem is present and the interpretation of the 


envelopes must be done with care, see the discussions in G raba rnik, Myllymaki , and 
Stoyan (2011) and in Baddeley, Diggle, Hardegen, Lawrence, Milne, and Nair 12014). 


Model evaluation is illustrated in the examples in Sections 3 and 4. 


3 Simulation experiments 


Realisations from the suggested models can be simulated sequentially using conditional 
distributions (|^ and (H-®. We recommend to use the scaled heterogeneity a{x) / max^g^y a{u) 
as the distribution for the first location, or alternatively, to condition to the first location 
Xi of data. Assume that (xi,..., Xk) are simulated. A simple algorithm for adding a point 
X to k = {xi,... ,Xk), or equivalently, simulating from the distribution having density 
fk+i{x\l^k), is to apply the accept-reject algorithm, see e.g. Ripley (1987, p. 61), which 
provides an upper bound for the conditional density. For the two models, the algorithm 
is as follows: 


• History-dependent rejection model: A point x following the conditional density 
a{x) K{xk,x) is proposed using the accept-reject method and the proposal is ac¬ 
cepted with the reweighting probability 7i(x, S{l^k, x)). 

• History-adapted model: The kernel width <pk = H{S{ltk)) is computed and pro¬ 
posals from the unnormalized transition density a{x) K^^{xk,x) are drawn using 
the accept-reject method. 


In the following simulation experiment we generate realisations of the new models with 
three parameter values in order to demonstrate the time evolution of the new processes 
and to see how data summaries capture their properties. We illustrate to what extent and 
how fast these processes can cover the target area, as well as whether the process starts 
to cluster in space. Each realisation consists of 100 points located in the unit square 
window. In this illustration, for the sake of simplicity, we assume that the target space is 
homogeneous, setting a{x) = 1, and hence the first point is drawn uniformly and the same 
starting point is used for all realisations. The first sampled point is xi = (0.22, 0.41). 
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3.1 History-dependent rejection model 


First, we demonstrate the history-dependent rejection model with self-interaction de- 
hned through the convex hull coverage, where the points outside the convex hull of the 
current point set are favoured according to the reweighting probability ([^. Second, we 
demonstrate the rejection model with recurrence self-interaction by using the reweighting 
probability (|^, which takes the number of generated points near the suggested point into 
account. 


3.1.1 Coverage self-interaction 


The purpose of this example is to illustrate self-interaction caused by the parameter 
p in the history-dependent rejection model with convex hull coverage (|^. We £x the 
parameter cx^ = 0.3 of the truncated Gaussian kernel ^ and vary the parameter p: 

• Model a, p = 1 (random walk without self-interaction). 

• Model 6, p = 0.1 (fast coverage), which accepts points inside the convex hull of 
previous points with low probability. 

• Model c, p = 0.5 (mild coverage), which accepts points inside the convex hull of 
previous points with mild probability. 


We simulate 19 realisations of the random walk model a, since it here represents a refer¬ 
ence model, and hve realisations of Model b and Model c. One of the simulated realisations 
of each model can be seen in Figure The polygons in the hgure illustrate the convex hull 
coverage related to the hrst 10 points of the process. One cannot detect much difference 
between the random walk model a and mild coverage model c, but the fast expansion of 
Model b can perhaps be seen: there are early (dark) points near the edges. 

However, point patterns only tell us about the spatial nature of the point process. Since 
we are mainly interested in the sequential (time order) aspect, we use four different 
functional summary statistics: ball union coverage (with radius 0.1), convex hull coverage, 
scanpath length and cumulative recurrence. (The two latter ones are explained below). 
The results related to these summaries can be found in Figure]^ The ball union coverage 
does not distinguish between the three models, but the convex hull coverage reveals that 
the coverage of Model b increases faster than for the other two models. Accordingly, 
Model b makes longer jumps on average than Model a or c, as can be seen from the 
scanpath length which measures the length of the sample path cumulatively. 

The recurrence function used in the reweighting probability (|^ calculates the numbers 
of points near the current point excluding the previous point, and the cumulative version 
sums all these numbers together. Now we see that Model b avoids the locations nearby 
other points when compared to the random walk. This may be due to the fact that the 
fast coverage process b penalizes slow coverage and hence increases the drift of points 
near the edges. 
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Model a 


Model b 


Model c 




Figure 2: Simulated patterns of the three history-dependent models with coverage self¬ 
interaction. The colour of the points denotes their order (from dark to light) and the 
first fixation is marked with a cross. The polygon indicates the convex hull of the first 
10 points. 






Figure 3: Ball union coverage with radius 0.1 (top left), convex hull coverage (top right), 
scanpath length (bottom left) and cumulative recurrence with radius 0.1 (bottom right) 
for the two models (hve realisations of each); The fast coverage model b is marked with 
black solid lines and the mild coverage model c is marked with red dashed lines. The grey 
area represents the envelopes estimated from 19 realisations of the random walk model a 
used as the reference model. 
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3.1.2 Recurrence self-interaction 


Here we illustrate self-interaction in the history-dependent rejection model with recur¬ 
rence We again fix the parameter of the truncated Gaussian kernel (|^ to 0.3 and 
vary the self-interaction parameter 6: 

• Model d, 6 = 0.5 (random walk without self-interaction). 

• Model e, 6 = 0.1 (low recurrence), which favours points in the non-visited areas 
rather than the points in the areas nearby the previous points. 

• Model f, 0 = 0.9 (high recurrence), which accepts points nearby the previous points 
with high probability. 


We again simulate 19 realisations of the random walk model d as well as five realisations 
of Model e and Model /. In Figure]^ we can see that, compared with the random walk 
model d, the realisation of the low recurrence model e indicates a tendency towards higher 
regularity, and the realisation of the high recurrence model / is clearly more clustered. 
Note also that Model e seems to cover the whole area quite fast compared with the other 
two models. 


Model d 



Model e Model f 




Figure 4: Simulated patterns of the three history-dependent models with recurrence self¬ 
interaction. The colour of the points denotes their order (from dark to light) and the 
first fixation is marked with a cross. The polygon indicates the convex hull of the first 
10 points. 

The results of the functional summary statistics are depicted in Figure The ball union 
coverage describes clustering of the points, hence for the high recurrence model e the ball 
union coverage curves locate lower than for the low recurrence model /, which covers the 
whole window with 100 points. The convex hull coverage curves reveal an effect similar 
to the ball union coverage: the low recurrence model e almost fills the whole unit square, 
whereas the high recurrence model / only fills around 60 % of the area. 

There is not much difference between the processes when comparing the scanpath lengths 
for the first 40 points, but after that the high recurrence model / makes shorter jumps 
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on average compared with the other two models. However, the cumulative recurrence 
function clearly reveals that the two models differ from the random walk model; the 
low recurrence model e avoids areas close to the previous included points, while the high 
recurrence model / favours areas near the previous included points. 






Figure 5: Ball union coverage with radius 0.1 (top left), convex hull coverage (top right), 
scanpath length (bottom left) and cumulative recurrence with radius 0.1 (bottom right) 
for the two models (five realisations of each): The low recurrence model e is marked with 
black solid line and the high recurrence model / is marked with red dashed line. The grey 
area represent the envelopes estimated from 19 realisations of the random walk model d. 


3.2 History-adapted model with convex hull coverage 
self-interaction 


In this example we fix the kernel width parameter r = 0.3 and pay attention to the effect 
of the parameter k of the history-adapted model (12)-(13). The parameter k controls the 
speed of decay as a function of coverage. We again define three history-adapted models: 


• Model g, K = Q (random walk, the kernel does not change in time). 

• Model h, K = 2 (mild clustering), which means that the process is allowed to take 
long jumps at the beginning, but eventually starts to cluster. 

• Model K = A (fast clustering), which starts to cluster rather quickly when the 
coverage increases. 
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The kernel function (12) here uses the convex hull coverage, which means that in 

is the area of the convex hull coverage generated by the points (a:i,..., Xk)- Now 
the history-adapted model works in such a way that at first the kernel width parameter 
r is dominating and the process can make long jumps, but when the area of the convex 
hull of points approaches the size of the window, the parameter k starts to affect and 
produces clustering. 


In Figures and it can be seen that the convex hull of the first 10 points is of the 
same size for all models: the speed of coverage seems to be similar for all processes at 
the beginning. The spatial structure of the mild clustering model h and the random walk 
model g are quite similar, but the points of the fast clustering model i are clearly more 
clustered than the points of the other two models, and there are only a few points in the 
upper half of the unit square. 


Model g 



Model h 



Model i 



Figure 6: Simulated patterns of the three history-adapted models. The colour of the 
points denotes their order (from dark to light) and the first fixation is marked with a 
cross. The polygon indicates the convex hull of the first 10 points. 

The functional summary statistics are plotted in Figure The fast clustering model 
i covers the area similarly to the other two models at the beginning, but after about 
50 points it starts to cluster, which can be seen as a decline of the ball union coverage 
summaries. The convex hull coverage does not reveal much difference between the models, 
and all the models are able to cover at least 60 % of the window. This is due to the wide 
kernel (r = 0.3) which allows the processes to make long jumps at the beginning. 

The summary statistic that shows the clearest difference between the three models is 
the scanpath length. While the jumps in the random walk model have time invariant 
transitions, the clustering models h and i start to take shorter jumps at some point, which 
is indicated by the decline in scanpath curves. In addition, the cumulative recurrence 
function shows that the fast clustering model i gathers points around the previous ones. 
To conclude, the effect of the decay parameter k seems to fasten the clustering as a 
function of coverage. 
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Figure 7: Ball union coverage with radius 0.1 (top left), convex hull coverage (top right), 
scanpath length (bottom left) and cumulative recurrence with radius 0.1 (bottom right) 
for the two models (hve realisations of each): The mild clustering model h is marked with 
black solid lines and the fast clustering model i is marked with red dashed lines. The grey 
area represents the envelopes estimated from 19 realisations of the random walk model g 
used as the reference model. 

4 A case study: Black Bow by Wassily Kandinsky 


We apply the developed modelling to experimental eye movement data related to arts in 
order to study self-interaction. The participants of the art experiment were inspecting six 
pictures of paintings and their eye movements were recorded. The stimulus picture was 
shown on the screen with a 1024 x 768 resolution and the eye movements were measured 
by the SMI iView X™Hi-Speed eye tracker with temporal resolution of 500 Hz. The 
distance between a participant’s head and the screen was about 85cm, and a forehead 
rest was used in order to prevent unintentional head movements. Each stimulus painting 
was shown for three minutes. The participants were also asked to describe the moods of 
the painting and their voice was recorded, but this information is not used here. 


We will focus on one painting used in the experiment, called Black Bow (1912) by Wassily 
Kadinsky shown in Figure (source of the painting: Diichting ( 1991[ )). We will £t 
four versions of the history-dependent model with recurrence self-interaction to the eye 
movement data of one subject. The goodness-of-ht of the model is checked using the four 
functional summary statistics mentioned earlier, and the best htting model is compared 
with the other subjects’ data in order to conclude whether the same model hts well for 
all participants. 
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4.1 Fitting the history-dependent rejection model with 
recurrence self-interaction 


We choose one subject of which eye movements are analysed and modelled here. Because 
of the long inspection period (three minutes), we decided to use only 100 hrst hxations of 
the sequence corresponding to a 35 second time-interval, as shown in Figure According 
to the two-stage model (see e.g. Locher et at, 1996) the overall impression of the scene 


is obtained during the hrst few hxations, and then the focus turns to the presumably 
interesting features. In addition, the gaze has a tendency to return to the interesting 
parts of the scene. Our aim is to hnd out whether we can hnd this sort of behaviour, i.e. 
if the process is of self-interacting type and if we can catch it with our rejection model. 


We hrst investigate the variation of the four functional summary statistics related to this 
particular data from Kadinsky’s painting. The ball union coverage with radius of 35 
pixels, convex hull coverage, scanpath length and cumulative recurrence (radius 50) of 
the 20 subjects of the experiment are presented in Figure [T^ as a dark solid curve for the 
subject under study, and as grey curves for the other participants. It can be noticed that 
the hrst 100 hxations do not cover the whole painting (the ball union covers around 30 % 
and the convex hull around 40 % of the target). It is typical of the eye movements that 
the edges of the painting are avoided and that is why the coverage hardly ever reaches 
the whole scene. 


Next, we estimate the heterogeneity term a{x) for the target painting. In this case, we 
utilize the empirical saliency map estimated as the intensity of hxation patterns of all 
the 20 subjects excluding the one under study (a total of 9366 hxations is used for the 
intensity estimation). Problems associated with the estimation of a{x) are considered 
in Discussion. For technical reasons a{x) is scaled to have values in [0,1]. The scaled 
saliency map together with the hxations of the subject under study can be seen in Figure 
This particular subject paid most attention to the areas with high intensity, but the 
gaze stayed in some low intensity areas also. 

In what follows, we ht four diherent rejection models consisting of three components: het¬ 
erogeneity (H), contextuality (C) and self-interaction (S). The self-interaction is assumed 
to be caused by the delayed recurrence function S'h( a^fc,a;), k = 2,... ,n — 1 dehned in 
([^, where the radius r = 50 pixels is used. Therefore, we condition by the hrst two 
hxations, Xi and X 2 - All four models are submodels of ([^, see Table 


Model 

Components 

Parameter values 

1 

H 

0 = |, large* 

2 

H, C 

6* = |, cr^ > 0 

3 

H, S 

0 < 6* < 1, cr^ large* 

4 

H, C, S 

0<9<l,a^>0 


Table 1: Components of the four rejection models. (*The value of cx^ should be chosen 
to be large enough such that the kernel is hat in the specihed window.) 


Model 1 includes heterogeneity and is a binomial process in a heterogeneous environment. 
When the Markovian kernel function is added (Model 2) the process is a random walk 
with Markovian property. Model 3 is a self-interacting process in a heterogeneous media 
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Empirical saliency map 



Figure 8; Empirical saliency map estimated from the fixations of all subjects excluding 
the one under study. The points indicate the hrst 100 fixation points of the particular 
subject under study. The hrst hxation is marked with a cross and the second with a star. 


without the contextuality effect and Model 4 contains both the Markov kernel and the 
self-interaction term. 


The log-likelihood function for Model 4 is now 


n—1 


n—1 


(16) 


k=2 k=2 

n—1 n—1 

+ log(l - 6') ^ l{Siiiltk,Xk+i)=0}i^k+l) + log(6') ^ 
k=2 k=2 

n—1 „ 

- log / «(w) ((1 - 0)l{S«(l^„n)=O} + 0 l{S«(l^„n)>l}) ^U. 

^=2 w 

The likelihood for Model 1 is just the hrst term on the right hand side of equation 
(16). For Model 2, the likelihood is obtained choosing 0 = 4 in (16) and for Model 3 
choosing cr^ to be large (i.e. large enough such that the kernel is hat in the specihed 
window). The parameters are estimated using the prohle-likelihood method with few 
iterative steps, see Section 2.3.1. for details. We used numerical integration for computing 
the normalization term of the log-likelihood (16), last row, and a grid of (60, 80,..., 400) 
and (0.05, 0.10,..., 0.95) for maximizing the likelihood. 


Model 1 includes only the empirical saliency map and we do not have to estimate any 
parameters. For Model 2 we get a = 180, and for Model 3 6 = 0.75. The parameter 
estimates for Model 4 are a = 180, 6 = 0.70. Note that for a random walk model we 
should have 6 = 0.50; hence there are recurrence features involved in this data. 


18 









4.1.1 Model comparisons 


Assessing the goodness-of-fit of the models is here done by estimating the four summary 
statistics mentioned earlier from the data and from 99 simulated realisations of the htted 
models. When simulating the model, we condition on the observed values of the hrst two 
hxations in order reduce variation right at the beginning of the process. One simulated 
realisation of each model can be seen in Figure and the summary statistics estimated 
from the data with pointwise envelopes estimated from the simulations in Figures [TO]- 
13 Note that the ball union and convex hull coverages are here presented with respect 


to the size of the window, hence they obtained values in [0,1]. 


Model 1 


Model 2 



Model 3 



Model 4 


r# 

L •*' V .. 4* 

. » •• 

- 

1 

, 


Figure 9: Simulated realisations of the four models htted to the eye movement data of 
the particular subject under study, overlaid with the empirical saliency map based on 
data from the other subjects. The hrst two points are hxed and marked with a cross and 
a star, respectively. 

For Model 1 all summary statistics show poor ht (Figure [l0|). Compared with the data this 
model covers the target area too fast, takes too long jumps according to the scanpath 
length, and goes to areas with too few points according to the cumulative recurrence 
function. As a conclusion, the heterogeneity component alone does not describe the data 
set well enough even though the spatial heterogeneity is followed rather well (see Figure 
upper left). 

Model 2 seems to perform slightly better than Model 1. The summaries estimated from 
the data set stay inside the simulated pointwise envelopes, except the ball union coverage 
and cumulative recurrence function after the hrst 70 points (Figure [TT|. These hndings 
indicate that data begin to cluster at the end of the inspection period, but the model 
does not carry that effect. 

Model 3 includes heterogeneity and self-interaction, but not contextuality, which is related 
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Figure 10; Ball union coverage with radius 35 (top left), convex hull coverage (top right), 
scanpath length (bottom left) and cumulative recurrence with radius 50 (bottom right) for 
the subject under study (solid line). Dashed lines represent pointwise envelopes estimated 
from 99 simulations of Model 1. 



0 20 40 60 80 100 0 20 40 60 80 100 

index index 


Figure 11; Ball union coverage with radius 35 (top left), convex hull coverage (top right), 
scanpath length (bottom left) and cumulative recurrence with radius 50 (bottom right) for 
the subject under study (solid line). Dashed lines represent pointwise envelopes estimated 
from 99 simulations of Model 2. 


20 













to the length of the jumps the process makes. As a result, this model seems to jump 
too much compared with data, since the estimated scanpath length summary is at odds 
with the simulated pointwise envelopes (Figure [I^. The marginal spatial structure looks 
slightly more clustered than Model 1 and Model 2 predict (Figure]^. 



index 


index 


Figure 12: Ball union coverage with radius 35 (top left), convex hull coverage (top right), 
scanpath length (bottom left) and cumulative recurrence with radius 50 (bottom right) for 
the subject under study (solid line). Dashed lines represent pointwise envelopes estimated 
from 99 simulations of Model 3. 

Model 4 includes all the three effects and seems to be in good agreement with data: all 
four summary statistics estimated for the subject under study stay within the simulated 
pointwise envelopes, see Figure [T^ It seems that this model is able to catch the nature 
of this eye movement process fairly well. The estimated parameter value 0 = 0.75 indi¬ 
cates that the locations nearby the previous points (excluding the most recent point) are 
favoured, which is a cause of spatial clustering. We conclude that the random walk model 
does not seem to be a good model for these data, but there is self-interaction due to the 
recurrence involved: the eye movement process seems to favour areas close to previous 
hxations. 


4.1.2 Population level comparison 

We have been able to describe the variation in the eye movement sequence of an individual 
by using the rejection model with recurrence-based weighting. In order to investigate the 
generality of the suggested models, we make comparisons at the population level using 
all 20 subjects. As can be seen in Figure [T^ the envelopes of Model 4 seem to cover 
the convex hull and ball union coverage curves of the subjects rather well. The scanpath 
length does not cover the curves as well: there are subjects whose gaze makes longer 
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Figure 13; Ball union coverage with radius 35 (top left), convex hull coverage (top right), 
scanpath length (bottom left) and cumulative recurrence with radius 50 (bottom right) 
for the subject under study (black solid line) and for all other subjects (grey solid lines). 
Dashed lines represent pointwise envelopes estimated from 99 simulations of Model 4. 


jumps at the beginning of the eye movement process than the fitted model predicts. 
Furthermore, the envelopes of the cumulative recurrence function cover almost all the 
curves, but there is one exceptional subject, whose fixations are strongly clustered after 
the 70th fixation. 


In order to further describe the variation in the data set related to self-interaction, we 
fitted Model 4 for the first 100 fixations of each subject separately. The estimates of the 
parameters a and 9 can be seen in Table The 90 % bootstrap confidence intervals for 
the parameter estimates are calculated from 20 realisations of the fitted model. When 
the parameter a is large, the model allows jumps over the observation window, and then 
the self-interaction parameter 9 dominates. Large 9 indicates strong spatial clustering. 
When a is small, one needs to take multiple jumps in order to cross the whole target 
window. 


In Figure 13, we have plotted the estimated summary statistics for all subjects together 
with the pointwise envelopes based on the model fitted to subject 5. For subjects 13, 
18 and 20, the scanpath length curves clearly exceed the envelopes of the Model 4 fitted 
for subject 5. For each of them, the fitted value of the parameter a in Model 4 is over 
280. This indicates that the process is allowed to make long jumps resulting in longer 
scanpaths. However, for these subjects the recurrence parameter 9 varies from 0.70 to 
0.80 and does not differ much from the recurrence parameter of subject 5. 


For subject 6, the cumulative recurrence curve is way above the envelopes, and the 
estimated recurrence parameter is 0.85 indicating strong clustering. The other subjects, 
whose cumulative recurrence curve is outside the envelopes, are 4, 12, 18 and 20, and for 
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Snbject id 

a 

90 % Cl for d 

9 

90 % Cl for 9 

1 

240 

200 

280 

0.75 

0.70 

0.85 

2 

180 

160 

200 

0.55 

0.45 

0.65 

3 

160 

140 

180 

0.65 

0.60 

0.75 

4 

160 

140 

180 

0.80 

0.75 

0.85 

5 

180 

140 

200 

0.70 

0.65 

0.85 

6 

180 

140 

200 

0.85 

0.80 

0.95 

7 

160 

140 

180 

0.70 

0.65 

0.80 

8 

240 

200 

280 

0.70 

0.65 

0.80 

9 

160 

140 

180 

0.65 

0.60 

0.75 

10 

160 

140 

180 

0.65 

0.60 

0.75 

11 

200 

180 

220 

0.70 

0.65 

0.80 

12 

260 

240 

320 

0.85 

0.80 

0.90 

13 

280 

220 

340 

0.70 

0.60 

0.85 

14 

220 

200 

240 

0.70 

0.60 

0.75 

15 

220 

180 

260 

0.55 

0.45 

0.60 

16 

200 

160 

220 

0.70 

0.65 

0.80 

17 

160 

140 

180 

0.70 

0.60 

0.80 

18 

340 

220 

420 

0.80 

0.75 

0.90 

19 

140 

120 

160 

0.70 

0.65 

0.80 

20 

280 

200 

340 

0.80 

0.70 

0.85 


Table 2; Estimated parameters a and 9 of Model 4 with their conhdence intervals for all 
snbjects. The snbject nnder closer stndy is nnmber 5 (bold). 


all of them the recnrrence parameter is over 0.80. As a conclnsion, there is some variation 
related to the clnstering effect of the points between these snbjects. However, for each 
snbject the estimated recnrrence parameter differs from 0.5 meaning that the random 
walk is not a snitable model. 


5 Discussion 


In this paper, we develop advanced data-analytical tools for extracting information from 
eye movement seqnences, which are needed in varions areas of application ntilizing eye 
tracking (see e.g. Rayner, 2009). Onr objective is to create simple bnt flexible dynamic 
stochastic models by employing mechanisms which nse the whole history of the seqnence 
in each gaze jnmp in order to captnre featnres of learning dnring the experiment. 


Heterogeneity of the scene, contextnality of snbseqnent fixations, and self-interaction 
of eye movements are elements that affect the eye movement process. We present a 
seqnential spatial point process approach which inclndes these three effects, the self¬ 
interaction being new in this context and is interpreted as a learning effect. This leads to 
what in probability theory is called self-interacting processes, which are generalizations of 
random walks in heterogeneons media. Althongh self-interacting random walks are well 
established in mathematics, physics and animal ecology, onr reasoning here is slightly 
different. We stndy how the process evolves at an early stage of an eye movement seqnence 
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whilst, e.g. in mathematics, the long term behaviour is of interest. Such processes are 
analytically difficult, even intractable, but their simulation is basically straightforward. 

After having constructed a new model, we need model fitting (estimation), evaluation 
of goodness-of-fit (model criticism) and simulation algorithms for various inferential pur¬ 
poses. Here, we suggest a likelihood approach for the new processes which is used in 
parameter estimation. It is not possible to obtain analytical results or use asymptotical 
reasoning. Instead, we compute the likelihood using simulation which allows us to make 
exact inference in the sense that it does not depend on the size of data and which in¬ 
cludes a boundary effect correction. In doing this, we enlarge the applicability of spatial 
statistics and the likelihood inference to a new area of applications. The processes can be 
used to make inference on the structure of data, including self-interaction, and to deduce 
uncertainties in conventional and new functional data summaries such as scanpath length 
and recurrence function. 


In this paper, we are interested in the dynamics of the eye movement process. The main 
question is whether there is self-interaction present in a given eye movement sequence 
and whether we can detect it using our new modelling. We focus on the beginning of an 
eye movement process, since, according to the two-stage model by Tocher and colleagues 
(see e.g. Tocher, [2006 ), the gist of the scene is established during the early fixations. 
Our history-dependent rejection model is, in fact, able to observe self-interaction in these 
particular data, although there is large within-subject variation. 


Our models utilize stochastic geometry in creating self-interaction caused either by cov¬ 
erage (how much of the scene is covered and how fast) or by recurrence (how much 
the process favours points nearby the previous points), both of which have justifications 
in eye movement literature. Functional summary statistics are needed for checking the 
goodness-of-fit of a fitted model, as well as for describing the structural components of 
the sequence. We use four summary statistics: convex hull coverage, ball union coverage, 
scanpath length, and cumulative recurrence. Several summary statistics are needed since 
none of these four was able to alone distinguish between the models in our simulation 
study. We found that the rejection model with convex hull coverage can be separated 
from the random walk by the scanpath length, whereas the rejection model with recur¬ 
rence can be distinguished from random walk with the coverage measures. The scanpath 
length and cumulative recurrence are needed with the history-adapted model for defining 
the speed of spatial clustering. 


We have illustrated two tractable process constructions for self-interaction, namely, history- 
based independent thinning and history-dependent transitions. These constructions are 
very different, and their use depends on the problem and the data set. The two devel¬ 
oped models are rather simple but can easily be extended. Other constructions are also 
possible, such as the heterogeneous mixture model, where the Markov kernel 
is replaced by 

p( affc) Ki{xk,x) + (1 -pi^k)) K2{xk,x). 

Here, Ki and K 2 are two kernel functions where the choice K 2 {xk,x) could be uniform 
in W or proportional to a{x), for example, and the mixture factor p(l^k) depends on 
the history of the sequence. Although simulation of such a model is straightforward, the 
associated inference is computationally demanding. 
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We have restricted our approach to sequential spatial point processes, mainly due to their 
tractability. However, this approach is a bridge to spatio-temporal models that would take 
fixation durations into account. For a separable spatio-temporal model, the spatial effect 
and time dynamics are multiplicative in the likelihood. A sequential spatial point process 
model can be used as a building block: if an order-dependent spatial model is available, 
the inclusion of time dynamics is straightforward, because inference on the ordered spatial 
aspect and fixation durations can be performed independently. If a preferred summary 
statistic contains information on the fixation durations, then a spatio-temporal model 
should be used instead of the sequential point process. This extension is a subject for a 
separate study. 


The estimation of the heterogeneity component a{x) is an issue not fully considered here. 
In the second order analysis of point patterns the first and second order components 
are not estimable from one observed point pattern without further information. Digglej 


et ah (2007) suggest two alternatives, which are the use of a parametric model for the 


intensity (heterogeneity) or, alternatively, the utilization of replications for the inten¬ 
sity estimation. In the sequential context the situation is similar. In our experiment 
sequences measured from several participants are available and are independent of the 
particular sequence under study. When using this information in the estimation of the 
heterogeneity component, the problem is that also these auxiliary control sequences are 
serially correlated leading to extra clustering at the sequence level. We assume that this 
effect is not as serious as in the intensity estimation from the case data only. When using 
auxiliary sequences, we have assumed that these sequences contain information which 
origins mainly from the target common to all participants and to the case under study 
and measures the wanted heterogeneity. 


An improvement would be to sample fixations from each of the auxiliary fixation sequences 
instead of using all the fixations as we did here. Merging these sampled fixations gives 
a point pattern which is then used in the estimation of a{x) using the kernel method. 
This will further reduce the effect of serial correlation. An alternative improvement 
is based on the case sequence under study by using the fitted model (containing both 
contextuality and self-interaction and a prefixed a{x)) to compute the inverse of the 
transition probability for each fixation. These weights can then be used in the estimation 
of a{x) by the weighted kernel method. The procedure can be iterated. The estimation of 
a{x) is discussed in Barthelme et al. (2013); Engbert et al. (2015). The reliable estimation 
of a{x) with open questions is a topic for a subsequent study. 


Another issue concerns the parameter estimation. Here, we suggest the profile likelihood 
or discretized coordinate descend algorithm for maximum likelihood using forward simu¬ 
lation. This early experimenting shows that by this method it is possible to separate the 
effect of self-interaction and present confidence intervals for parameter estimates and con¬ 
fidence envelopes for chosen summary statistics. We know that approximative inference, 
being computationally much faster, is also a possibility and would be very important 
in a methodological toolbox. The extensive experimenting with approximative inference 
related to this is a future task. 
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